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| j Abstract 

We investigate equilibrium properties of small world networks, in which both connectivity and spin vari- 
i | ables are dynamic, using replicated transfer matrices within the replica symmetric approximation. Popula- 
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£h . tion dynamics techniques allow us to examine order parameters of our system at total equilibrium, probing 
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both spin- and graph-statistics. Of these, interestingly, the degree distribution is found to acquire a Poisson- 
like form (both within and outside the ordered phase). Comparison with Glauber simulations confirms our 
results satisfactorily. 
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I. INTRODUCTION 



Small worlds are systems characterized by a local neighborhood (given by short-range bonds) with 
a sparse set of long-range connections per spin. This simple architectural effect has been shown to 
bring about remarkable cooperative and synchronization phenomena. The term 'small-world' has 
been coined by the now famous experiment by the Harvard social psychologist Stanley Milgram [jj: 
in 1967, as part of his research on the network of acquaintances in the United States, he took a 
number of letters and handed them over to people totally unrelated with the addressees, and, with 
the instructions to pass them over to someone they think might know the addressee. This process 
was repeated until the letters finally arrived to their destination. Stanley Milgram then estimated 
the average path length from the two randomly chosen individuals which turned out to be a mere six. 
This experiment revealed that although social networks are very sparse, in reality any two pair of 
nodes can be topologically very close. In fact, numerical studies of other types of real networks (e.g. 
citation, linguistic, disease spreading, etc) show that the small-world effect is a common architecture 
among real network structures and brings about optimal information processing. The question then 
arises, how do networks spontaneously evolve from (almost) random configurations into particular 
structures such as small-world ones? And which underlying process drives the distribution of the 



long-range short-cuts within the nodes? The above questions fal 



under a particularly active area 



of research, namely the evolution of networks (see e.g. DD , or for a recent review). Since real 
networks (be it biological, social, economic or otherwise) hardly ever maintain a static architecture 
this problem of predicting network structure has important applications. In this paper we attempt 
to formulate and describe the thermodynamics of the problem from an analytic point of view. This 
carries the obvious set of advantages and disadvantages: while resulting in robust and exact results, 
it will be amenable to a set of (perhaps not fully realistic) assumptions. To be precise, we examine 
a coupled system on a small-world architecture in which both nodes and connections are mobile. 
However, the two dynamic processes occur on distinct timescales; connections are assumed to evolve 
slowly enough such that, at each of their update steps, spins have effectively reached equilibrium. 
This will allow us to avoid solving the explicit dynamical relations and instead turn directly to the 
thermodynamics. Our starting point is the free energy per connection degree of freedom. We couple 
the two dynamic processes of the spins and the connections by constructing two Hamiltonians: a 
typical Ising one describing the energy of the spins and a Hamiltonian of the connections, constructed 



to reward network configurations minimizing the free energy of the spins. This choice allows us to 
proceed analytically while retaining a sufficient amount of realism. The result is a replica theory 
where the replica dimension represents the ratio between the two temperatures (of the spins vs 
connection processes). 

Our paper is organised as follows: In the following section we introduce our model and the 
pair of energy functions describing the thermodynamics of the spins and the graph variables. In 
section UTTl we write the total free energy of the system as an extremisation problem in terms of the 
typical finite-connectivity order parameter function. We then proceed to define the observables of 
our system, of which there are here two kinds, probing spin (section II V A J) and graph (section HVBjl 
organization statistics respectively. The single pure state approximation (section E| allows us to 
deal with the resulting replicated transfer matrices following the diagonalization process of j^, 0]. 
We first derive in section IV Al numerically tractable forms for our set of order functions which are 
to be solved via population dynamics. Observables such as magnetization, average connectivity, or 
degree distribution then follow easily, see section IV Bl We perform a bifurcation analysis and plot 
phase diagrams showing the transition lines between ordered and paramagnetic phases in section 
IV CI We find that, perhaps contrary to initial expectations, the resulting degree distributions are 
close to, or exactly, Poisson. Comparison with numerical simulations shows good agreement given 
the complexity of these experiments requiring adiabatic (practically infinitely long) timescales. 



II. MODEL DESCRIPTION 



We study a system of A" Ising spins cr = (<7i, . . . , a^) with <jj G { — 1, 1}, arranged on a "small- 
world" structure. We represent this by a one-dimensional lattice with uniform nearest-neighbor 
interactions of strength J and with randomly-chosen sparse short-cuts of strength Jy G {—J, J} 
that can connect distant pairs of spins (i, j). For every i ^ j we assign a variable Cjj denoting 
whether a connection exists (cjj = 1) or not (cjj = 0), with cu = 0. In the absence of short-cuts the 
average path length is N/4 while in the combined system the scaling is bounded above by log(N). 
This significant reduction in the path length is commonly termed the "small- world" effect 7]. For 
static architectures, in which the link and bond matrices {cij, J^} are taken as quenched random 
variables, frustration effects are known to induce spin-glass phases p|. Our model aims to examine 
thermodynamic properties of the above spin systems under the freedom of allowing the connectivity 
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and bond matrices {cy, Jy} to evolve in time in search of the state that best promotes order within 
the system. To be precise, on short timescales the links and bonds can be seen as static variables 
with respect to which the spins equilibrate, while on longer timescales cy and J y - explore their 
configuration space. The measure of this latter process is related to the ordering within the spin 
system on the instantaneous state of the graph. Thus the spins and the graph architecture on which 
they live are dynamically interwoven. It is quite natural that the architecture dynamics depends 
on the entire system state (including the spins) rather than just the architecture itself (as with e.g. 
preferential attachment jsj]). Links and connectivities are taken here to evolve on identical timescales, 
although generalizing this to more involved scenarios is also possible. On the timescale in which spins 
reach thermal equilibrium our combined system is described by the "fast" Hamiltonian 

H f (a, c, J) = -J ^2 a i a i+i ~ Z~2 a '- l 'J r 'i a J ( 1 ) 

i i<j 

(where we take periodic boundary conditions on the chain). Spins equilibrate with respect to (JIJ at 
an inverse temperature 1//?/, and their behavior is described by the partition function 

Z f (c 1 J) = Y,e-' 3fHfi(T ' C ' J] (2) 
cr 

On timescales sufficiently long to guarantee that spins have reached equilibrium, links and bonds 
are not static, but evolve dynamically and we will take their stationary state to be described by the 
"slow" Hamiltonian 

H s (c,J) = — -j- log Zf(c, J) + V(c, J) (3) 

Pf 



V(c,J) = jz2 Ci i log (ir) +logcosh(K p )-/i- p ^ 



(4) 

\ C I J 

Kj 

This choice energetically favors those configurations of {cy, Jy} that minimize the free energy of the 
spins. The role of the chemical potential V(c, J) is twofold: firstly, it aims to preserve the overall 
nature of the small-world system; it guarantees that for N — > oo the average number of connections 
per spin is a finite number. Secondly, it allows us to tune the relative concentration of { — J, J} bonds 
in the system (as we will see in section IIVB1 the former is controlled by the variable c whereas the 
latter by K p ). The connectivity- and bond-variables {cy, J^} equilibrate with respect to this slow 
Hamiltonian at inverse temperature f3 s , leading to a total partition function 

Z S = Y, e-^ c ' J ) = £)[2>(c, J)]^' e-^ c - J ) (5) 
c,J c,J 

4 



This partition function, by construction, contains n = /3 s //3f replicas of the fast system. In general, 
the ratio of inverse temperatures n can take any value (integer or otherwise) so that analytic con- 
tinuation in the replica dimension depends solely on our choice of temperature values. The limit 
n — > corresponds to temperatures T s —>■ oo in which the partition sum (jJJJ) is dominated by the 
entropy of the slow system. In contrast, T s — > favors those architectures {c^, Jy} that increase 
order among the spin variables for a given number of links. Note that this is a general optimization 
criterion which does not enforce a priori any particular structure on the links but allows the links 
to arrange themselves. In fact, the graph statistics become interesting observables, which we can 
measure, rather than enforced constraints. Our order parameters follow from the slow free energy 
per spin 

fs — — h m - — log Z s (6) 
and derivatives of this generating function. 

III. THE FREE ENERGY 

To calculate the slow partition function (jSJ) we first take the trace over the connectivity and bond 
variables {cy, Jy} 

e^' JoS ' <Tl - cri + 1 Y[ [l + — (e /3fJOV<Tj ) J (7) 



<J 1 ...(T N i<j 



up to irrelevant multiplicative constants. We denote cr, = (aj, . . . , erf) where <Tj • crj = J2 a a f a< j an d 
defined the abbreviation (f(J))j = [2 cosh(i^ p )] _1 [e Kp /(J) + e~ Kp f(—J)]. We are interested in the 
case where c (the chemical potential for bonds) is finite, and hence c/N — > in the limit iV — > oo, 
so that the above product can alternatively be seen as a product over exponentials (up to terms of 
0(N~ 2 )). We thus encounter the typical nested exponential form of finite connectivity problems. 
To achieve site factorization it is convenient to introduce into our expressions the order parameter 
function P. 

i 

via appropriately defined delta functions, which is a probability distribution over replicated spins. 
In the limit iV — > oo we can now evaluate the free energy (JHJ) via steepest descent and express it as 
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an extremization problem in the space of probability distributions -P(cr), namely: 

k ^ s crcr' ^ tT!...a N i > 

where Tcr ,cr' [-P] represent the transfer matrix elements 



Tcr.er't-P] = ex P 



(9) 



and -P(cr) is to be evaluated from the fixed-point equation 



Tr (Q\<j] T N \P\\ 

P{°) = k {t n [p]) Qcr,cr' [t] = 8<r tT 6v0> (10) 



For more details on the derivation of the above expressions we refer the reader to |y| where the special 
case of the n — > limit was studied. 

Finding solutions of (fTUj) amounts to diagonalising the transfer matrix T of dimensionality 2 n x 2 n . 
This problem has been solved in jsj. Here we will not be concerned in the entire spectrum of 
eigenvalues, as the limit N —>■ oo ensures that only the largest eigenvalue Ao will provide a non- 
vanishing contribution to the free energy. The left- and right- eigenvectors associated with this 
eigenvalue follow from the equations 



J2 T <T,<T>iP]U(<T')=\ U(*) (11) 

cr' 

Y, V (<r') T <T>M p ] =A F(cr) (12) 



cr 



These eigenvectors are unique up to the usual arbitrary multiplicative factor, and non-negative 
[5, We note that we need both left- and right-eigenvectors since the transfer matrix T[P] is 
non-symmetric. The order function -P(cr) (jHj) is manifestly normalized. Due to our scaling freedom 
for the eigenvectors we can always choose them so that J2cr ^( cr ) = J2cr V(. CT ) = 1- The physics of 
our system is given by the normalized distributions -P(cr), V(er), U(cr) which are to be found by self- 
consistently solving equations (JOHEJ) (in fact, U(cr) and V(er) turn out to represent the distributions 
of cavity spins with a chain bond rather than a graph bond removed [ill]). 
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IV. OBSERVABLES 



A. Spin system observables 

We are interested in probing organizational properties of our system both within the spin variables 
and the connectivity ones. For the spin system, we define the canonical observables; the magnetiza- 
tion and the overlap order parameter as moments of the probability distribution (jHJ), namely 

m a = ^P(cr) a a (13) 

q aP = ^P(<T)aV (14) 

(T 

In the above and henceforth, the quantities P(tr), V(cr), U(a) are given by their saddle-point values. 

It is well known that infinite dimensional systems, such as small world lattices, with frozen bonds of 
random signs, will have a spin-glass ground state at low temperatures for certain values of the control 



parameters 
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13, 



14|. This spin glass ordering is intimately linked to frustration within the system; 
the inability of spins to find energetically optimal configurations. By allowing the architecture some 
limited degree of freedom, we expect that the system will be able to optimize its state somewhat 
better. Probing the degree of frustration within the system as the slow temperature is varied is 
therefore an interesting problem. The frustration is normally defined as the fraction of closed loops 
from sites ii — > i 2 — — *■ *i where the product J ili2 . . . J ikil is negative. Unfortunately, to 
measure this directly in our system where bonds are mobile would require us to be able to measure 
correlations over long length scales within the system (in fact scaling like the average loop length 
~ log(iV)), which is technically difficult. To try and finesse this problem, in 3, Qj- the fraction of 
misaligned spins was calculated, i.e. the fraction of spins that did not point in the direction of their 
local field. Due to the mobility of the connections in our system we expect that thermal equilibrium 
states within the ordered phases will be steered towards configurations where spin alignment with 
their local fields is optimal. The result of this structural organization can be measured by the 
following quantity = J"° dh P(l, h) + dh P{— 1, h) which gives the fraction of misaligned spins 
and is defined in terms of the joint spin-field distribution P(a, h) = liniAr^oo i 'Yl li {5a,a i 5[h — hi(<r)]) g 
where (• ■ -) s denotes thermal averages over the slow process (x) s = Z^ 1 j e~@ sHs ( c > and 
hi(a) = 52 - CijJijdj + Jo(<Jj + i + (7j_i) denotes the local field at site i. However, at e.g. very low 
temperatures, one expects the spins to align to their local fields whether they are in a spin glass phase 
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or not. Thus to try and get a different measure to probe the frustration in the system we consider 
the fraction of bonds in the graph, which are not energetically optimized by the spin configuration: 

^ = ( Q (- a i a i+i J o)) s + -tfj ^2(CijQ(-^jJij))s (15) 

i i<j 

This is also not an absolute measure of frustration, but in the low temperature spin glass phase i/j 
will be non-zero, as opposed to a low temperature ferromagnet where we would have ip = 0. The 
calculation of either ^ or is similar to the calculation of the free-energy, with a specific observable 
(i.e. matrix in the transfer matrix notation) at one or two sites. We find 

V> = D 1 J2 V(<r)e(-aio[J )T trir ,[P\U(* r ) + D 2 J^ P( ( T)P( ( r')(Q(-a 1 a' 1 J)e^ J(T - (T ')j (16) 
fro" era" 

where D\ and D 2 are normalization constants to give the fraction of sites, i.e. D\ = 

V(<r)7W[PM«r') and D 2 = £ ff<r P{a) P{<t>) {e^™'). 



B. Connectivity system observables 

Let us now inspect organisational phenomena within the graph. We first identify the roles played 
by the control parameters c and K p that appear in the chemical potential (@J). This can be done by 
adding suitable generating terms into the Hamiltonian Q and monitoring their impact on For 
instance, if one transforms H s — >• H s + Yl,i<j c ij then taking derivatives ^|a=o translates to 

^^E^')— c E p ( CT ) p ( T ) <^ J<TT > J (it) 

ij <TT 

Now c represents the average number of connections per spin in our system. It depends on the 
replica dimension n via the scalar spin product and it reduces to c = c in the limit n —>■ 0. In the 
limit c — > oo (scaling J as J/c to keep the local fields in the graph hf(cr) = Ylij c ij^ij a j °f ^(1)) 
we again recover c = c to leading order as found in [15j. Similarly to the above, one also finds that 
taking H s — > H s + ^J2i<j c ijJ-ij produces the average bond strength on the graph. As well as being 
interested in the above average connectivity and bond strength at total equilibrium, we would also 
like to investigate the connecti vity structure in more detail. To make contact with a variety of recent 
work on complex networks 0, ^| we define the degree distribution for our system 

i 
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Following a calculation similar to that of the free energy in section UTTl one easily finds that 

/c\h 
™ e* kk V(<t)U(<t') exp[cJ2(P(r) e^^j + {3 f J cr ■ <r'] (19) 
n crcr' r 

The above observables are all expressed in terms of the trio of distributions P{cr), V(cr) and U(cr), 
taken at the saddle-point of the free energy (JHJ). To proceed with a numerical evaluation of the 
observables one now needs to specify a form for these. 

V. REPLICA SYMMETRY AND TRANSFER-MATRIX DIAGONALISATION 

To solve the self-consistent equation (JTUJl one is required to make certain assumptions. Firstly 
with regards to the form of the order function P(cr) and the eigenvectors. They represent different 
distributions over replicated spins (for any We will consider the simplest possible scenario 

in which permutation of spins within different replica groups a = 1, . . .,n leave the order function 
invariant (replica symmetry). This is equivalent to assuming the existence of a single pure state. 

For any natural n G N + it is relatively straightforward to express these distributions, as their 
support is a finite discrete set. However, for the more general case of n G K one has to make an 
analytic continuation which leads to more complicated expressions. 

For any natural n we can impose replica symmetry by writing for any arbitrary distribution X{&) 

n 

X{p) = J2 X ( £ ) S[2£-n;J2 a c,} (20) 

e=o a 

where normalisation of X(ar) requires Yle=o ^Ww) = !• O n the other hand, for any real n 

X(*) = fdz x(z) f[ 6 T (21) 

where now normalization requires J dzx(z) = 1. The above ansatze hold for any distribution X(er), 
and in particular as X G {P, U, V} we define the natural n ansatze in terms of {V,U, V} and for real 
n in terms of {p, u, v} respectively. 

A. Self-consistent order function equations 



Our self-consistent equations for P, U, V (jl()H12)) can now be transformed into relations between 
the field distributions {V,U,V} or {p, u, v}. Let us start with the natural n versions. It is conve- 
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nient to begin by working out an identity for the replica symmetric form of the general expression 
X(er)F(cr ■ r). We insert the replica symmetric ansatz (|2Uj) for X, use the gauge transformation 
°a - ► o" q t q and introduce the representation of unity 1 = J2k=o — ra ! J2 a T °<\ wmcri results in 

n n 

x{*)f{ct • -o = e E x ^w k - E r J E ^ 2£ - n; E ^] F (E ^) ( 22 ) 

£T f=0 fc=0 a <7 a a 

We now define the set of replica indices S = {a £ {1, . . . , ra} : r a = 1} and its complement 5 = {a6 
{1, . . . , n} : r a = — 1} which allows us to write r a a a = J2 a &s — ^2aes Isolating these last 
two summations via the unities 1 = Ylki=o <^[2&i — Sags an d 1 = Sfc 2 ~=o ^[2^2 + k — n; J2aes a °\ 
and using the general identity ^ CTi ap 8[2q — p; J2a=i a <*\ = (q) j we obtain 

n fc n—k 

Y,x{*)f{* . r) = E E *0 ^ 2fc - n ' E r J 

CT £,fe=0 fci=0 fc 2 =0 a 

x 5[£ + k + k 2 - h - n; 0] F(2(k 1 + k 2 ) - n) ( * ^ ( " ~ ^ (23) 



\ &2 / 

Using the above identity (and very similar manipulations) we can write our self-consistent equations 

as 

U{£) = X^in) exp[cA P (e, J)] A v (£, J ) (24) 

n j n-j 

V(£) = Aq \n) E E E V ^ e cAp(j ' J)+p f Jol2(kl+k2) ~ n] 5[£ +j + k 2 -k 1 -n; 0] (25) 

j=0 ki=0 k 2 =0 

V{1) = U ^ V ^ (26) 

The largest eigenvalue Ao(ra) follows from the above by utilizing the normalization condition 
Q/^CO = 1- We have introduced the convenient shorthand, 

n l n-l /A / — A 

A x (£, J) = E E E X ® ^ + k + 1 ~ 3 - "5 0] ( ) ( n k ) (e^ +k ^)j (27) 

for X £ {P, U, V} and X £ {V,U, V} respectively 

Let us now turn to the more general case of ra £ IR. Firstly, a Taylor expansion of the transfer 
matrix elements © into a series of exponentials and insertion of the replica symmetric ansatz (j21jl 
leads to 

Tv^P] = ePf^-V ( e ^E„^) e ( 28 ) 
10 



with (•••)$ representing averages over the measure 




and where we introduced the functions 



i<k i {J[} 



A(J,x) = atanh(tanh(/3/ J) tanh(x)) (30) 
B(J,x) = ^log[4cosh[/3/J + x]cosh[/3/J-x]] (31) 

The first of the above functions can be identified as a 'message' (or effective field) passed during 
)elief propagation while the latter is related to the free energy shifts which occur during an iteration 
I4I \v\ . For n — > and within replica symmetry this second term does not contribute although 
in the more general case of n > it will play an important role. Following the belief propagation 
picture, one can also relate (j2*9Jl to a weighted measure over the messages coming from the long range 
bonds. Performing the spin summations in (|11H12|) using the ansatz (|2Tjl and requiring the resulting 
expression to have the eigenvector form leads to 

A (n) u(x\n) = [ dx' u{x'\n) ^ J 3 ? U nB ^^5[x -9- A{J , x')]) (32) 

A (n) v(y\n) = J dy' v(y'\n) (e"^' +e )% - A(J , y' + 0)]) g (33) 

so that the largest eigenvalue follows from the above by simple integration. To close the above 
equations we also need to derive an expression for the function p(h). The starting point for this is 
equation (fTUj). Rewriting the traces in terms of the eigenvectors and substituting our ansatz (j2*Tj) 
results in 



f dx dy u{x) v(y) < 


cosh(h) 


n 

\ S[h 


- (x + y)] 


2 cosh(a;) cosh(y) 


J dx dy u{x) v(y) < 


cosh(x+y) 


r 


2 cosh(x) cosh(j/) 



P\ h ) = r — ^ ( 34 ) 



The coupled set of equations ()32|33|34|) are to be solved self-consistently. This has a clear inter- 
pretation in terms of message-passing algorithms: p(h) gives the distribution of messages passed 
along long-range short-cuts, whereas, u(x) and v(y) that of messages passed along the chain (from 
;he left- and right- neighbors). This inspires a solution using a population dynamics methodology 



171 ]: the main difference is that one is now required to weight the averages of the field distributions 
by an n-dependent factor. In practice, expressions of the form <fi{x') = J dx 4>(x) w(x) S[x' — g(x)} 
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(for some arbitrary probability density 4>(x), weight w(x) and updating function g(x)) are solved by 
sampling values of x from the density (f>{x) and updating x' — ^ g(x) with weight w(x). To interpret 
this weighting term, one can write w(x) = \_w(x)\ +p, where, ^(^J is the integer part of w(x), and 
p is the fractional part. At each iteration step we replace ^(^J °f the population members with x' 
and a further member with probability p. 



B. Observables within replica symmetry 

To calculate the magnetization (j!3)l and spin glass $11}) order parameter within our system, we 
substitute the replica symmetric ansatz for P{<r) (|2()|) into their definitions, which together with a 
minor rearrangement gives for integer n 

™=±Mr_ i y n ±M n - 1 ) (35) 

while for real n the expressions become: 

m = J dh p{h) tanh(/i) (37) 
q = dhp(h) tanh 2 (/i) (38) 



So given V or p from the self-consistent equations, for either integer or real n, we may evaluate these 
order parameters. In figure^we plot the magnetization for two different values of n, and compare our 



results against simulation experiments. More details on the simulations are given in section IV D[ but 
we note here that these experiments are particularly time consuming due to the coupled dynamics, 
so that, only modest system sizes are allowed for reasonable CPU cost. Within these constrains we 
feel that the agreement is reasonable. 

Evaluating the fraction of energetically non-optimal bonds ip is slightly more involved. We use 
the replica symmetric representation of the transfer matrix which for integer n reads 

n 

T a(T , [P] = $( 2k - n ;J2 a ^ e p f Jo(T - (T ' +cAp{k ' J) (39) 

fc=0 a 
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1.25 



0.75 



m 0.5 



0.25 




-0.25 



T f /J 

FIG. 1: We plot the magnetization m as a function of temperature Tf for n = 1 and n = 5. The solid lines 
are the theoretical predictions while the dotted lines are a guide for the eye joining the markers with error 
bars which come from simulations. We have c = 2 and Jo = Jij = 1 Vi, j. The simulations were done via 
Monte Carlo Glauber dynamics (see text for details) onjV = 200 spins. Despite the small system size they 
seem to be in reasonable agreement with the theory. 

Then, after some combinatorial work, ip is found to be given by 
^ = D, E E * PfM2k - j - 1] f? ~ ! )(* ~ 1 V? ~ I ) \ V(i)U(jy Ap{l ' J) + V(j)U(i)e cA ^ 



=1 j=0 k=0 

n n—l i—1 



n n—i i—L / i \ / 

+ fl 2 2 P v:5:5:nopw(": 1 1 )( 

i=l j=0 fc=0 ^ ' ^ 

i=l j=l fc=0 

77. — 1 Tt — 1 % 



i — 1 J \ k J \j — k J 
n — 1\ /i — 1\ (n — i 



k j \ j — k 
n — 1\ /« — 1 



^ } J(n+2{2k-i-j)) 



i - 1 
n—l 



n — i 



k J \j — k — 1 
n — i — 1 



3 /3 / J(n+2(2fc+l-i-j)) 



i=0 j=0 fc=0 

For real n the replica symmetric transfer matrix is given by, 



(3 f J(n+2(2k-i-j)) 



leading to 

i/> = Di / dxdyd^ u(a;)u(y)M(0) 



To-^[P] = e^ Jo<r -°" y d6M(6) e ^ 6 ^^ 

e'Pf 10 2 cosh (x + 6-y) 



(40) 



(41) 



+D 2 / d/iid/i 2 p(hi)p(h 2 ) e"^ J 



[e^2 cosh(a; + 6» + /3/ J ) + e~^2 cosh(x + 6- P/J )] 
Q(J)2 cosh(/ii + h 2 ) + Q(-J)2 cosh(/ii - h 2 ) 



n-l 



[e h2 2 cosh(/i + p f J) + e-^2 cosh(/i - /3/J)] n_1 
13 



(42) 



Tf/J 



FIG. 2: We plot the fraction of misaligned bonds ip against the re-scaled temperature Tf/J for c = 2 and 



Jo = 1. The solid lines are for n = 1 while the dotted lines are for n = 3. The upper pair of lines are for 
r = 0.8 while the lower pair are for the ferromagnet, r = 1. We see that in the ordered phase, increasing n 
allows the system to optimize the bonds energetically. 

where we define 



We plot ip for a few sets of parameters in figure El As the fast temperature goes to oo we have 
ip — > 1, i.e. exactly half of the bonds at any point in time are energetically optimal, so in the high 
temperature phase the ordering is non-existent (as one would expect). We also see that increasing 
n leads to better levels of optimization, again what we would expect but it is possible to quantify it 
here. Decreasing r and hence increasing the disorder, makes it harder for the spins to energetically 
optimize themselves, although at low temperatures, due to the condensation phenomena in the bonds 
(see below), the magnetization will increase to 1. With all spins aligned, the fraction of energetically 
non-optimal bonds becomes exactly the fraction of bonds with J < 0. We see this in figure El as 
Tf — > 0, ip — ► 1 — r. This would not be the case for the pure spin-glass, r = 0.5. 

We now turn our attention to the graph observables. We first focus on the average connectivity 
which is expressed as 




r = 



2 cosh(K p J) 



(43) 




(44) 
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FIG. 3: We plot the average number of bonds c against temperature for J = Jq = 1 and c = 2 in the 
ferromagnet. The higher line with the first order transition is for n = 5 while the lower line is for n = 1. 
The dotted lines are the theoretical predictions while the solid line is a guide for the eye linking the error 
bars which are measurements from simulation experiments. The agreement is reasonable although, as in fig. 
^ we find that the sharp transition is smeared out due to the small system size for our simulations (N = 
200 spins). 

for integer n and 



If figure El we plot the average number of bonds c. At low temperatures (the specific temperature 
depends on other parameters) the average connectivity increases sharply. This is due to ordering 
within the spin system leading to increased energetic gain by adding connections. Higher values of n, 
for a given Tf, means a lower value of T s and hence the connectivity variables will be governed more 
strictly by the free energy of the fast system, which is minimized by high connectivity configurations. 

We also looked at the full degree distribution, which is given up to normalization constants for 
integer n by 



c = c dhi dh,2 p(hi) p(h 2 ) {cosh(/3/ J) + tanh(/ii) tanh(/i 2 ) sinh(/3/ J)} 



n 



(45) 



for real n. 




(46) 
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FIG. 4: We plot the probability that a given node has degree k, Pk against k for n = 0.5, r = 0.6, 
Tf = J = Jq = 1 and c = 2 within the spin-glass regime where other observable values are m = and 
q ~ 0.581. The true degree distribution is given by crosses, for comparison we have also given the Poisson 
distribution with the same value of c with circles. Although there are differences between the two for this 
set of parameters, the difference is very small. 



and for real n by, 

EW ~ if 



n 

£<k 



dh e dJ£p(hi)Q(J e 
[2 cosh(/^)]" 



dxdyu(x)v(y) 
[4 cosh(a;) cosh(?/)] T 



^nY li B(J l ,h l )+nB(J ,x) 



(47) 



x 



2 cosh 



+ A{J , x)+Yl A ( J *> h * 



A typical example of this degree distribution is given in figure |U What is particularly interesting is 
that although the degree distribution is in principle free to take on any form it keeps very close to 
that of the Poisson degree distribution with mean c. In fact in the paramagnetic phase, we know 
that P{cr) = 2~ n and thus we find c exactly from (|17p without invoking replica symmetry, namely 
cpm — ccosh n (f3J), which is independent of r since cosh is an even function. Thus the average 
degree is independent of the bond disorder (in this model) in the paramagnetic phase. Here, the 
degree distribution also scales linearly with c. By using the fact that in the paramagnetic phase 
we also have U(cr) = V(cr) = 1~ n we can also see that S(fc) = e~ CpM Cp M /k\, i.e. the degree 
distribution is exactly Poisson. We also find exact results in the fully ferromagnetic phase where 
P(er) = U{cr) = V{cr) = EL-Wi- There c FM = c(e^ nJ )j and S(jfe) = e-~ CFM c k FM /k\. In both these 
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cases the degree distribution is exactly Poisson. We can understand this, since in both phases there 
is no energetic gain in having any particular c^- = 1, since it will not affect the spin distribution 
(they are either all set to be aligned, or fully random) and thus the degree distribution will be the 
maximum entropy one, i.e. Poisson. It is also clear that there are a range of ordered states between 
the two extremes above, and that we cannot say anything further analytically about the degree 
distribution there. However, we may evaluate our order parameter equations numerically, and we 
find that although the degree distribution is not Poisson, it is very close, as shown in figure EJ It was 
not obvious that this should be the case, and indeed, the increased critical temperature for a scale 
free degree distribution would have suggested that this could be optimal, since it increases ordering, 
but it transpires that this is not the case here. 

C. Phase diagrams 

Having derived the main equations from which our observables follow, we can now proceed to 
the evaluation of the transition lines in our phase diagram numerically and via a bifurcation anal- 
ysis. Firstly, we see that the state p(x) = u(x) = v(x) = 5(x) always solves equations (j32|33|34j) . 
giving m = q = 0, for all temperatures. We can therefore associate this state with the high- 
temperature (paramagnetic) solution. To examine continuous bifurcations away from this solution 
we assume that close to the transition the fields are small and that the paramagnetic ^-distributions 
evolve to either distributions of small, non-zero mean (in leading order) marking the paramag- 
netic/ferromagnetic transition or to distributions of small, non-zero variance (again in leading order) 
marking the paramagnetic/spin-glass transition. With these considerations in mind we define the 
moments h l = J dhp(h) h e = 0(e e ) for some < e C 1 (and similarly for x l = J dxu(x)x e and 
y e = J dyv(y)y e ). We assume that there is no first order transition. Then, expanding equations 
()32l33l34j) for small values of fields and using h = x + y and h 2 = x 2 + y 2 which follows from (}3"4*j) we 
arrive at paramagnetic/ferromagnetic and paramagnetic/spin-glass transition lines: 

P -> F : 1 = c (sinh^/J) cosh™" 1 ^ J)) J e 2f3 f J ° (48) 
P -> SG : 1 = c (sinh 2 (/3/J) cosh n " 2 (/3 / J)) J cosh(2/3/ J ) (49) 

These equations reduce to those found in [6] in the limit n — > 0, recovering the small-world bifur- 
cations. The correspondence is exact if we identify the paramagnetic mean connectivity here with 
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FIG. 5: We plot the phase diagrams for c = 2 and J = Jo = 1. The left figure is for n = 0.1, where 
the solid lines are given by the bifurcation conditions (|48J49[) while the markers come from solving the 
order parameter equations numerically and the dotted line linking markers is a guide to the eye. The right 
figure is for n = 2 and all lines are linking markers which come from solving the order parameter equations 
numerically. The P— >F and P— >SG transitions are here first-order. For larger values of n we see that the 
links are better able to align to increase order; firstly, the transition temperature from the paramagnetic 
phase is higher and secondly, the size of the spin glass phase is significantly smaller. 



c. It also reduces to those of [la] for J = and Q(J') = S(J f — J), where the Hopfield model on a 
dynamic random graph was studied, if in the latter only a single pattern is stored (in this scenario 
the Hopfield model becomes equivalent to a ferromagnet with a different gauge). It is well known for 
these models [n| that as n increases, the transitions are increasingly likely to be first order. Thus 
to produce phase diagrams of the system, as well as looking at the bifurcation lines given by the 
above we also solved the full equations numerically. The results are shown in figure |3] where we see 
that increasing n decreases the size of the spin glass phase, which we expect is due to the increased 
cooperativity. 
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D. Simulations 



In order to check the validity of our theoretical work we performed numerical simulations of this 
model. To do this we need to introduce a dynamical process on both the spins and the graph which 
will converge to an equilibrium distribution described by their respective Hamiltonians. One way to 
do this is via Glauber dynamics |18| . the dynamics then automatically obey detailed balance. The 
transition rates between a given state and another state with a single "spin" flip (where we take spin 
in the broader sense to include the binary variables {cy} and { J^} as well as the more familiar {a - ;}) 
is determined by half the energy difference (or local field) between the two states. Defining general 
spin flip operators via F^$(cn, Cu, . . . , Cnn) = $(cn, Cu, . . . , —Qj, . . . , Cnn) and similarly for 
the Glauber rates can be written as 



W[F? jC ,c] = ^1-tanh 



2 N 2 



(50) 

W[Fp,J] = I|l-tanh[j^ p -^log(e- 2 ^ c - J ^)]} (51) 

where the angular brackets denote averages over the fast process for the given realization of the graph 
and bonds. 

The nature of the coupled dynamics means that for each change to the graph (the slow dynamics) 
one must re-equilibrate the spins, measure the averages as required in the above equation and sub- 
sequently change the graph configuration again. Thus the computation effort required to equilibrate 
the slow system is very large compared to simulations on a given, fixed, graph. In particular, for 
strongly disordered graphs, where changing a single bond is expected to seriously alter the free-energy 
surface, it is very difficult to obtain reasonable statistics. Instead, we have focused our efforts on the 
simpler case of purely ferromagnetic bonds. This means that after changing a given bond, the new 
equilibrium distribution is expected to be very close to the old one and also the equilibration times 
will be shorter in general. We have performed simulations on systems with N = 200, and in figures 
^ and H3 we compare the results with our theoretical predictions. Due to the small system size we 
must expect that there are both persistent errors due to the relatively small system size, smearing 
of all phase transitions and large error bars on any given measurement. Bearing this all in mind we 
feel that the results, particularly for the average connectivity, clearly support the theory. 
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VI. CONCLUSIONS 



The study of complex networks has recently become a very popular field due to their ubiquity in 
nature, technology and social interaction, where these fields are meant in a broad sense. While the 
statistical structure characterizing real world networks (path lengths, degree distributions,. . .) and 
models that recreate these properties have been extensively studied from experimental measurements 
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We 



idv J 

3- 



on real world systems, through numerical simulations and theory, understanding the be 
networked systems based on local rules (dynamics) is still a relatively unexplored area 
have presented a solvable model that examines a spin system on a small world graph with which 
we have probed cooperative behavior of the entire system (both of the graph and the spins). To 
overcome the theoretical challenge of systems evolving on disparate timescales we have focussed on 
the adiabatic limit; the graph evolves infinitely slowly relative to the spin variables. This allows us 
to treat the model using the well developed thermodynamics of replica theory, rather than having 
to treat the dynamics explicitly. The advantage of this approach is twofold. Firstly the results are 
exact in the thermodynamic limit in the region where replica symmetry is stable. Although we have 
not examined replica symmetry breaking experience suggests that this would only occur for n < 1, 
at low temperatures (high values of (3f) and for some critical amount of disorder in the bonds { Jij}- 
The second benefit of this approach is related to the relative simplicity of our present approach. 
We do not specify in advance the dynamics of the graph, but instead only describe it through its 
equilibrium energy function. Thus the resulting graph structure becomes an observable itself, rather 
than an object which is fixed from the start. Indeed, naive intuition may suggest that the optimal 
structure could have been scale free, so that ordering in the spins would have occurred at a higher 
temperature. It turns out that this was not the case, apparently due to entropic reasons. 
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